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calculation.  In  the  calculations  these  fluctuations  arise  from  small  numerical  inaccuracies. 
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the  incident  shock  wave  leading  to  non-uniform  reflection,  thermal  conduction  to  trie 
walls,  and  interactions  with  boundary  layers. 
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WEAK  AND  STRONG  IGNITION 

I.  NUMERICAL  SIMULATIONS  OF  SHOCK  TUBE  EXPERIMENTS 


I.  Introduction 

Anomalously  low  and  somewhac  erratic  ignition  temperatures  in  mix¬ 
tures  of  hydrogen  and  oxygen  were  first  reported  by  a  number  of  authors 
doing  shock  tube  or  adiabatic  compression  experiments  [1,2, 3, 4],  Later 
the  work  of  Strehlow  and  Soloukhin  and  their  coworkers  [5,6,7,8,9,10] 
described  variations  in  the  ignition  behavior  of  mixtures  ignited  behind 
reflected  shocks.  Strehlow  catalogued  three  types  of  behavior  behind  the 
reflected  shock:  (1)  acceleration  of  a  reflected  shock  when  a  pressure 
wave  generated  by  released  chemical  energy  reaches  the  reflected  shock; 
(2)  an  accelerating  wave  headed  by  a  shock  behind  the  reflected  wave; 
and  (3)  a  fast  release  of  energy  leading  to  a  detonation  behind  the 
reflected  wave.  Voevodsky  and  Soloukhin  [10] ,  who  concentrated  on  varia¬ 
tions  in  ignition  due  to  chemical  kinetics,  defined  two  ignition  regimes 
separated  by  an  extended  second  explosion  limit.  "Strong"  or  "sharp" 
ignition  was  said  to  occur  when  the  results  of  schlieren  streak  photo¬ 
graphs  showed  that  a  single  ignition  locus  at  the  reflecting  wall  was 
sufficient  to  form  a  detonation  front.  "Mild"  or  "weak"  ignition 
occured  at  lower  temperatures  where  the  reaction  was  initiated  at  many 
loci  which  finally  merge  together  to  form  a  uniform  front. 

The  work  of  Meyer  and  Oppenheim  [11]  showed  that  mild  ignition 
started  in  distinct  centers  in  the  shock  tube,  usually  in  stagnant 
corners.  They  then  showed  that  the  chemical  induction  time  can  be  very 
sensitive  to  temperature  variations  which  occur  due  to  nonuniformities 
behind  shock  waves.  Subsequently  they  developed  a  "coherence"  theory 
for  ignition  which  gave  criteria  for  strong  ignition  based  on  a  suffi¬ 
ciently  high  rate  of  chemical  energy  release  in  a  large  enough  radius 
Manuaeript  submitted  September  15, 1981. 
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[12,13].  Properties  of  these  "ignition  kernals"  or  "exothermic  centers" 
were  described  by  Zajac  and  Oppenheim  [14],  van  Tiggelen  [15]  and  Borisov 

[16] .  Borisov  described  the  origins  of  the  exothermic  centers  as  re¬ 
sulting  from  temperature  fluctuations,  fluctuations  is  activated  mole¬ 
cules  or  radicals,  pressure  waves,  or  catalytic  generation  of  radicals 
on  suspended  particles.  He  pointed  out  that  it  is  virtually  impossible 
to  eliminate  turbulence  and  pressure  waves  due  to  non-uniformities  in 
temperature  and  concentration.  These  fluctuations  cause  the  induction 
times  of  reactive  mixtures  to  vary  in  time  and  space. 

In  this  paper  we  use  data  from  shock  tube  experiments  by  Cohen  and 
Larsen  [17]  combined  with  numerical  simulations  to  study  weak  and  strong 
ignition.  As  Borisov  suggested,  we  find  that  the  weak  ignition  system 
is  extremely  sensitive  to  small  acoustic  perturbations.  In  these  cases, 
the  detailed  one-dimensional  simulations  described  below  show  fluctuation 
behavior  quite  similar  to  that  seen  in  the  multi-dimensional  experiments. 
When  the  systems  are  not  as  sensitive,  as  in  the  strong  ignition  case, 
we  find  using  a  previously  validated  set  of  rate  constants,  that  the  si¬ 
mulations  and  the  experiment  are  in  excellent  quantitative  agreement. 

II.  The  Experiments 

The  calculations  described  in  this  paper  are  based  on  two  of  a 
series  of  reflected  shock  tube  experiments  performed  by  Cohen  and  Larsen 

[17] .  The  experimental  technique  consisted  of  simultaneous  measurements 
of  pressure,  UV  light  transmission  (either  emission  or  a  combination  of 
emission  plus  absorption) ,  and  recording  through  time  resolved  schlieren 
photography  of  the  shock  reflection-reaction  wave  formation  process. 

The  shock  tube  and  technique  for  taking  time-resolved  schlieren  photo¬ 
graphs  have  been  described  by  Strehlow  and  Cohen  [15].  Pressure  was 
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measured  by  a  fast  response  piezoelectric  transducer  (3  usee  risetime) 
flush  mounted  at  the  center  of  the  end  wall  of  the  shock  tube.  The 
transducer  signal  was  recorded  by  means  of  an  amplifier  (Kistler  Model 
565)  and  an  oscilloscope.  Another  pressure  transducer  was  located 
198  mm  upstream  from  the  back  wall,  primarily  as  a  trigger  for  oscillo¬ 
scopes  and  counters. 

For  the  light  transmission  measurements,  a  monochromator  (Schoeffel 
Model  QPM30)  was  used  to  isolate  a  region  centered  at  3100°  A  (theore¬ 
tical  half  band  width  =  50°  A) .  A  photomultiplier  (1P28)  monitored  the 
transmitted  UV  light  intensity.  The  effective  absorption  path  is  re¬ 
stricted  to  a  region  1  mm  from  the  back  wall  and  1  mm  wide  by  slits  lo¬ 
cated  on  the  monochromator  and  the  shock  tube  window.  The  signal  was 
recorded  by  means  of  a  cathode  follower  and  oscilloscope.  The  risetime 
of  the  electronic  circuits  was  less  than  1  ysec.  Figure  1  contains  a 
schematic  diagram  of  the  experimental  arrangement. 

The  typical  experimental  cases  will  be  studied  in  detail  using  the 
model  described  below.  One  set  of  initial  conditions  is  in  the  strong 
ignition  limit  and  the  data  is  compared  directly  to  the  results  of  a 
simulation.  The  other  experimental  conditions  are  in  the  "mild"  or 
"weak"  ignition  limit.  Because  the  timescales  for  this  latter  case 
are  relatively  long,  the  cost  of  a  simulation  of  these  exact  para¬ 
meters  is  correspondingly  high.  Thus  we  have  simulated  a  problem  which 
is  still  in  the  weak  ignition  limit,  but  the  temperature  and  pressure 
have  been  raised  slightly.  The  results  of  the  simulation  are  then  used 
to  guide  the  interpretation  of  the  experiment.  Figures  2  and  3  show  the 
schlieren  photographs  for  the  experimental  parameters  described  in  Tables 
1  and  2. 
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Table  I 


Parameters  for  the  Strong  Ignition  Experiment 


H2:02:Ar/2:l:7 

Undisturbed 

Incident 

Reflected 

Temperature 

298  K 

621  K. 

1036  K 

Pressure 

0.066  atm 

0.362  atm 

1.3  atm 

Fluid  Velocity 

4 

4.76x10  cm/ sec 

Shock  Velocity 

7.54xl04  cm/s 

4 

4 . 5x10  cm/  s 

Mach  Number 

2.165 

Table  II 

Parameters  for  the  Weak  Ignition  Experiment 

H2:02:Ar/8:2:90 


Undisturbed 

Incident 

Reflected 

Temperature 

298  K 

601  K 

972  K 

Pressure 

0.233  atm 

1.07  atm 

3.4  atm 

Fluid  Velocity 

3.74xl04  cm/s 

Shock  Velocity 

6.59xl04  cm/s 

4.36xl04 

Mach  Number 


1.98 


III.  The  Numerical  Model 

The  conservation  equations  for  mass,  momentum  and  energy  [18]  solved  in 
the  calculations  performed  below  may  be  written 


dp 

at 


V 


•pv 


(1) 


V  *  n.V.  -  V  •  n.v  +  Q 
-  1 


L.n. 

J  J 


(2) 


=  -  V  •  (pvv)  -  VP 


d*  =  _  v  .  Ev  _  7.  Pv  -  •  £  W 

at  -  ~ 

where  the  heat  flux,  Q,  is  defined  as 


Q  =  -  aVT  +  Z  njhjXj 

j 


(5) 


The  quantities  p,  pv,  &  and  P  are  the  total  mass,  momentum,  energy  density, 
and  pressure,  respectively.  Here  v  is  the  fluid  velocity,  and  k  is  Boltmann's 
constant.  The  {n^}  and  {V^}  are  the  number  density  and  the  diffusion  velo¬ 
cities  of  the  individual  chemical  species.  The  quantity  X  is  the  thermal 
conductivity  coefficient  of  the  gas  mixture  at  specified  {n_. }  and  tempera¬ 
ture,  T.  The  (Qj)  and  refer  to  chemical  production  and  loss  processes 

for  species  j .  The  last  term  in  Equation  (5)  represents  the  local  change 
in  energy  due  to  molecular  diffusion  and  chemical  reactions  which  must  be 
added  to  the  fluid  dynamic  energy  density.  The  quantities  (tu)  are  the 
temperature  dependent  enthalpies  for  each  species.  We  note  that  viscosity 
and  thermal  diffusion  have  been  omitted  from  Equations  (l)-(5)  since  they 
have  a  negligible  effect  on  the  solutions  presented  below. 
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The  reactive  shock  model  used  to  perform  the  calculations  described 
below  solves  these  equations  in  one  spatial  dimension  for  the  generalized 
coordinate,  r.  Though  r  may  represent  position  in  either  cartesian,  cy¬ 
lindrical  spherical  or  some  generalized  coordinates,  the  particular  applica¬ 
tion  here  is  in  cartesian  geometry.  The  ideal  gas  law  is  assumed  so  that 
P  -  NkgT  .  (6) 

where  N  is  the  total  number  density, 

N  *  L,  n.  .  (7) 

j  3 

Thus,  we  may  write  the  internal  energy  per  unit  volume  as 

e  =  L  n.h  -  P  -  H  -  P  (8) 

j  J  n 

j 

The  internal  energy  is  related  to  the  pressure  and  total  energy  by 


P  «  ]  2  _  1  n 

e  -  — r  =  G  -  -r  ov  -2j  h  .n.  s  E  -  —  pv“ 

Y— 1  2  j  o]  J  2 

where  y  is  the  ratio  of  specific  heats,  {hQ^ }  is  the  set  of  heats  of 

formation  at  0°K  of  the  species  j ,  and  E  is  the  total  energy  minus  the 

total  heat  of  formation.  The  diffusion  velocities  {V.}  are  found  by 

-3 

inverting  the  following  matrix  equations  [18,  19]: 


(9) 


“A 


-  v. >  =  E  w4U  a  -  v.)  , 


i  k  n2d  -*  ~j 

jk 


jk  ^-k  -4c' 


(10) 


subject  to  the  constraint 


P.  V. 
J  ~j 


0 


The  source  terms  are  defined  as 

S  £  7  (n./N)  -  (p.  h  -  n./N) 
J  “  J  3  3 


VP 

P 


(ID 


\ 

4-/ 


k 


n  n,  _ 

±-  lVVW 


SD 


jk 


VT 

T 


.(12) 


Throughout  this  paper  physical  quantities  are  generally  give  in  cgs  K  units 
although  pressure  is  also  given  in  atmosphere. 
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The  convective  and  diffusive  transport  terms  in  Equations  (l)-(4) 
are  solved  separately  and  then  coupled  together  by  timestep  splitting 
with  asymptotics  as  described  by  Oran  and  Boris  [ isJ •  The  convection 
terms  are  solved  using  one  variant  of  the  Flux-Corrected  Transport  (FCT) 
method  [20,21].  This  is  a  conservative,  monatonic  algorithm  with  fourth 
order  accuracy  which  does  not  require  artificial  viscosity  to  stabilize 
shocks.  The  ordinary  differential  equations  describing  the  chemical 


kinetics  are  represented  by  the  and  the  L^n  terms  in  Equation  (2). 
These  are  solved  using  VSAIM,  a  version  of  the  CHEMEQ  algorithm  [22,23], 


which  has  been  specially  optimized  for  use  in  dynamic  models.  For  most 


of  the  calculations  presented  below  the  timescales  under  consideration 
are  short  and  the  diffusive  transport  processes,  thermal  conduction  and 


molecular  diffusion,  have  negligible  effect.  Thus  the  calculations 


performed  omit  these  processes  unless  otherwise  indicated. 


The  detailed  simulations  of  the  shock  tube  experiments  require  that 
we  generalize  the  adaptive  gridding  method  used  previously  [24]  to  in¬ 
clude  two  finely  gridded  regions.  The  first  region  moves  with  the  shock 
front,  surrounding  first  the  incident  and  then  reflected  shock  front.  The 
second  region  stays  at  the  reflecting  wall  until  the  reaction  wave  has 
developed  and  then  moves  with  the  reaction  wave.  Each  region  may  have  a 
different  minimum  computational  cell  size.  When  the  shock  and  reaction 
waves  merge,  the  detonation  front  then  maintains  the  cell  size  characteris¬ 
tic  of  the  most  finely  zoned  region.  For  many  of  the  calculations  presented 
below,  computational  cell  sizes  varied  nearly  two  orders  of  magnitude  from 
the  most  finely  zoned  region  around  the  reaction  wave  to  the  coarsest 
zoning  in  the  calculation. 


Figure  4  shows  che  geometry  and  general  configuration  of  the  calcula¬ 
tions  used  to  simulate  the  experiments.  The  calculation  is  done  in  Cartesian 
geometry  and  initialized  with  an  incident  shock  moving  from  right  to  left. 

The  left  hand  boundary  is  a  reflecting  wall,  and  the  right  hand  boundary 
condition  ensures  that  there  is  a  constant  inflow  of  material  moving  at 
the  incident  shock  velocity.  We  thus  minimized  computational  costs  by  not 
modelling  the  full  problem  which  includes  a  driver  section  of  the  shock  tube. 
This  is  valid  because  throughout  the  calculated  time  span  the  reflected 
shock  wave  remains  well  separated  from  the  contact  discontinuity  formed  when 
the  diaphragm  initially  burst. 

The  chemical  kinetics  rate  scheme  used  consists  of  about  fifty  chemical 
rates  relating  the  species  H^,  0^,  0,  H,  OH,  HO^,  H^O^,  an<*  H2^'  ^aS 
been  tested  extensively  against  experimental  data  [25,  26,  27,  28]  and  shown 
to  give  excellent  results.  Recent  calculations  have  shown  that  it  also  gives 
values  of  the  flame  velocity  which  compare  well  with  experimental  measure¬ 
ments.  The  reaction  rate  scheme  is  give  in  Table  3.  Heats  of  formation 
and  enthalpies  have  been  taken  from  the  JANAF  tables  [29]. 
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Table  III.  H,-0,  Elementary  Reactive  Mechanism 


tt 

•y-4 

B 

AT  exp 

(-C/T)(a) 

Reaction 

A(h) 

B 

c(b> 

Refere.ce 

H  4  HO  -  0  +  H 

1.40(-L4) 

L  .00 

3.50(403) 

[331 

*■  2 

3.00(-14) 

1.00 

4.48(403) 

[33] 

H  4  HO-  "*•  H,  4  0, 

4 . 20(-ll) 

0.00 

3.50(402) 

[331 

2  +-  2  l 

9.10(-11) 

0.00 

2.91(404) 

[33] 

H  +  HO,  *  HO  +  HO 

4.20(-10) 

0.00 

9.50(+02) 

[331 

2.00(-ll) 

0.00 

2.02C+04) 

[33] 

H  4  H02  *  0  4  H20 

8 . 30(-ll) 

0.00 

5.00(+02) 

[34] 

1.75(-12) 

0.45 

2.84(+04) 

k  -  k,/  K 
r  f  c 

H  4  HO,  *  HO,  4  H2 

2 . 80(-12) 

0.00 

1.90(403) 

[33] 

1.20(-12) 

0.00 

9.40(403) 

[33] 

H  +  H,0,  -  HO  4  HO 

5.28(-10) 

0.00 

4.50(403) 

[33] 

2  2  ^  2 

3.99(-lO) 

0.00 

4.05(404) 

k  -  k,/K 
r  f  c 

HO  4  H,  "•  H  4  H,0 

1.83(-15) 

1.30 

1.84(403) 

[35] 

1.79(-14) 

1.20 

9.61(403) 

[35] 

HO  +  HO  ~  h2  +  02 

1.09(-13) 

2.82(-U) 

0.26 

0.00 

1.47(404) 

2.42(404) 

k  -  k-/K 
f  [36] 

HO  4  HO  —  0  4  H.O 

1.00(-L6) 

1.30 

0.00(400) 

[35] 

**■  2 

3,20(-15) 

1.16 

8.77(403) 

k  -  k./K 
r  f  c 

HO  4  HO,  -  H,0  +  0, 

8. 30(-ll) 

0.00 

5.03(402) 

[37] 

2.3S(-10) 

0.17 

3.69(404) 

k  -  k./K 
r  f  c 

HO  4  H„0  *  HO,  +  H, 

1.70(-il) 

0.00 

9.10(402) 

[33] 

Z  «-  2  2 

4 . 70(-ll) 

0.00 

1.65(404) 

[33] 

HO  4  0,  -  HO,  4  0, 

1.60(-12) 

0.00 

9.56(402) 

[34] 

3-^2  2 

6.69 (—14) 

0.33 

2.04(404) 

kr  -  y*c 

HO  4  H,  -  HO  +  H,0 

1.20(-12) 

0.00 

9.41(403) 

[36] 

2  ■*-  2 

1.33(-14) 

0.43 

3.62(404) 

k  -  k./K 
r  f  c 

HO,  4  HO,  -  H,02  +  0, 

3.00(-ll) 
1.57 (-09) 

0.00 

-0.38 

5.00(402) 

2.20(404) 

[34] 

k  -  ke/K 
r  f  o 

Table  III.  (continued) 
H.,-0.,  Elementary  Reactive  Mechanism 


ki  -  ATB  exp  (-C/T)(a) 

Reaction  Reference 


A(b) 

8 

c(b) 

0  -  HO  -  K  +  On 

2 . 72 (-12) 

3. 70(-10) 

0.28 

0.00 

-8.10(+01) 

8.45(403) 

kf 

-  k  /:; 

r  c 

133] 

0  4  HO,  -  HO  4  0 

8 . 32(-ll) 

0.00 

5.03(402) 

[37] 

2.20(-ll) 

0.18 

2.82(404) 

k 

r 

»  k'K 
f  c 

0  4  -  H20  +  02 

1.40(-12) 

0.00 

2.12(403) 

[34] 

5. 70(-14) 

0.52 

4.48(404) 

k 

r 

-  k./K 
f  c 

0  +  H.,0.,  -  HO  +  HO., 

1.40(-12) 

0.00 

2.13(403) 

[34] 

2.07 (-15) 

0.64 

8.23(403) 

k 

r 

-  VKc 

H  +  H+  M—  H-  4  M 

1.80 (-30) 

-1.00 

0.00(400) 

[33] 

^  fa 

3. 70(-10) 

0.00 

4. 83 (-04) 

[33] 

H  +  HO  +  M  -  H,0  4  M 

6.20(-26) 

-2.00 

0.00(400) 

[331 

5. 80 (-09) 

0.00 

5.29(404) 

[33] 

H  +  0,  +  M  -  HO  +  M 

4. 14(-33) 

0.00 

-5.00(402) 

[33] 

*-  *“  fa 

3. 50 (-09) 

0.00 

2.30(404) 

[33] 

HO  4  HO  +  M  ;  H202  +  M 

2.50(-33) 

0.00 

-2.55(403) 

[331 

2. 00 (-07) 

0.00 

2.29(404) 

[33] 

0  +  H4M  —  H0  4M 

8.28(-29) 

-1.00 

0.00(400) 

[38] 

2.33(-10) 

0.21 

5.10(404) 

k 

r 

-  k./K 
f  c 

0  +  HO  4  M  -  HO.,  +  M 

2.80(-31) 

0.00 

0.00(400) 

[38] 

1.10(-04) 

-0.43 

3.22(404) 

k 

r 

*  k./K 
f  c 

O+O+M-O,  +M 

5 . 20(-35) 

0.00 

-9.00(402) 

[33] 

3. 00 (-06) 

-1.00 

5. 94(404) 

[33] 

(a)  ‘Umolecular  reaction  rate  constants  are  given  in  units  of  cm^/(molecule  sec). 

Teraolecular  reaction  rate  constants  are  given  in  units  of  cm  /(molecule”  sec). 


(b)  Exponentials  to  the  base  10  are  given  in  parenthesis;  i.e.,  1.00(-10)  » 
1.00  X  10"10. 


IV.  Simulation  of  the  Strong  Ignition  Experiment 

The  initial  calculations  performed  tested  the  individual  chemical 
and  hydrodynamics  algorithms  separately  [30],  First  we  performed  a  purely 
chemical  kinetic  calculation  to  determine  the  induction  time  of  a  mixture 
of  H^sO^sAr/Z: 1: 7  at  the  reflected  shock  temperature  and  pressure.  Figure  5 
shows  a  series  of  a  constant  volume  adiabatic  calculations  performed  for 
various  initial  temperatures  using  the  CHEMEQ  algorithm.  The  induction 
time  for  the  purpose  of  this  paper  is  defined  as  that  time  at  which  the 
temperature  increases  just  20  K  from  its  initial  value.  As  can  be  seen 
from  the  figure,  the  induction  time  at  1034K  is  109+  1  ysec,  in  good  agree¬ 
ment  with  the  times  as  measured  by  the  shock  tube  experiment.  Similarly, 
tests  of  the  hydrodynamics  algorithm  alone  also  give  shock  velocities  in 
agreement  with  experiment.  Thus  we  know  that  the  individual  chemical  and 
hydrodynamic  parts  of  the  model  are  correct.  The  next  step  is  then  to  use 
the  model  to  study  the  interaction  of  these  processes  in  the  reactive  flow. 

One  problem  with  numerical  simulation  even  in  one  dimension  is 
interpreting  the  large  quantity  of  output  data.  At  each  timestep  we  cal¬ 
culate  mass,  momentum,  energy,  temperature,  and  a  number  of  species  densities. 
In  order  to  compare  the  simulation  most  effectively  with  the  experiment,  we 
have  chosen  to  follow  the  time-dependent  behavior  of  two  quantities:  the 
maximum  pressure  in  the  system  modelled  and  the  fluid  velocity  at  the  re¬ 
flecting  wall.  The  maximum  pressure  at  each  timestep  is  constant  until  the 
incident  shock  is  reflected;  then  it  jumps  to  a  higher  value.  It  increases 
again  when  chemical  energy  is  released  and  the  reaction  wave  starts,  and 

reaches  a  maximum  at  the  time  when  the  reaction  wave  and  reflected 
shock  wave  merge.  The  fluid  velocity  behind  the  reflected  shock  in  the 
ideal  case  is  zero.  In  the  simulation  it  is  zero  until  the  incident  shock 
hits  the  wall.  At  this  point  it  oscillates  wildly  to  maximum  values  of  nearly 
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Fig.  5  —  Calculations  of  chemical  induction  time  as  a  function 
of  temperature  for  the  strong  ignition  case 
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1x10^  cm/sec  and  then  settles  quickly  to  a  very  small  value  until  the 
reaction  wave  starts  and  it  jumps  up  again.  Then  the  fluid  velocity  near 
the  wall  quickly  decays  again,  and  rises  to  a  substantial  value  when  the 
wave  generated  from  the  merging  of  the  shock  and  reaction  wave  reaches 
the  wall. 

Figure  6  shows  a  tracing  of  the  maximum  system  pressure  and  the  wall 
velocity  from  the  computer-generated  output.  There  we  see  that  the  time 
calculated  for  the  formation  of  the  reaction  wave  is  in  good  agreement 
with  the  experiment.  There  is  apparent  disagreement  between  the  calcu¬ 
lation  and  experiment  in  the  individual  times  calculated  for  merging 
of  the  shock  and  reaction  wave  and  the  reflection  of  this  wave  from  the 
wall,  although  the  sum  of  these  numbers  is  very  close  to  that  determined 
experimentally.  In  fact,  there  is  some  ambiguity  (perhaps  five  micro¬ 
seconds)  in  determining  these  quantities  from  the  original  schlieren 
photograph . 

Since  the  overall  behavior  of  the  simulation  appears  to  agree  well 
with  the  experiment,  we  look  now  at  the  results  in  more  detail.  Figure  7 
shows  the  temperature  and  pressure  of  the  reaction-wave  before  it  merges 
with  the  reflected  shock  front.  It  looks  like  a  detonation  propagating 
into  the  shocked  mixture.  In  fact  it  has  become  a  detonation  at  about 
35  Usee  after  the  reaction  wave  starts,  as  seen  in  Figure  8.  Figure  8 
shows  the  positions  of  the  shock  front,  reaction  wave,  and  merged  wave 
and  the  contact  discontinuity  formed  as  the  waves  merge  as  a  function 
of  time.  The  reaction  wave  starts  out  slowly  and  then  accelerates  to 
a  detonation  travelling  at  a  velocity  which  is  2%  less  than  the  Chapman- 
Jouet  velocity  and  about  15%  greater  than  the  experimental  velocity. 

After  it  merges  with  the  reflected  shock  front,  it  decelerates  relative 
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TIME  (mmc)-  —  ► 

Fig.  6  —  Comparison  of  calculation  and  experiment  for  the  string  ignition  case. 
Top  panel:  schlieren  photograph  with  relative  times  marked.  Bottom  panel:  cal¬ 
culations  of  mavirnmn  pressure  in  the  system  and  fluid  velocity  at  1  mm  from 
the  wall  as  a  function  of  time  with  relative  times  marked. 
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Transmitted 
Detonation ' 


Fig.  8  —  Calculated  position  of  the  reflected  shock  front  reaction  wave, 
transmitted  detonation  and  contact  discontinuity  as  a  function  of  time 
for  the  strong  ignition  case 


to  laboratory  coordinates  due  to  the  incoming  incident  shock  flow.  It 
moves  at  a  velocity  which  is  about  13%  less  than  the  Chapman- Jou get  velo¬ 
city  and  about  6%  greater  than  the  observed  experimentally.  The  contact 
discontinuity  formed  when  the  waves  merge  moves  forward  more  slowly,  at 
the  fluid  velocity.  Figure  9  shows  the  temperature  profiles  after  the 
waves  have  merged  for  calculations  using  two  different  grid  resolutions, 
indicating  that  the  calculations  have  converged  numerically. 

Finally,  Figure  10  shows  the  calculated  number  density  of  OH  at  a 
distance  of  1  mm  from  the  reflecting  wall.  In  the  experiment  shown  in 
Figure  2,  OH  was  first  detected  by  the  absorption  experiment  85  usee 

after  shock  reflection.  Thus  we  see  from  the  graph  that  the  number  den- 

3  -3 

sity  at  this  time  is  about  2x10  cm  .  The  observed  emission  delay  time 

for  OH  was  126  usee;  at  this  time  in  the  calculation  the  OH  number  density 
17  -3 

is  approximately  10  cm 

We  conclude  then  that  for  this  case  studied,  the  agreement  between  the 
experiment  and  the  simulation  is  remarkably  good.  This  is  especially  gra¬ 
tifying  since  it  is  a  direct  test  of  both  the  individual  components  of 
the  calculation  and  of  the  procedure  used  to  couple  the  chemical  kinetics 
and  the  convective  transport.  The  calculation  is  so  good  that  one  is 
tempted  to  use  it  to  calibrate  the  experimental  diagnostics. 

V.  Simulation  of  a  Weak  Ignition  Case 

Historically  the  weak  ignition  regime  in  reflected  shock  experiments 
is  one  in  which  ignition  is  unpredictable.  A  neat,  clean,  distinct  reaction 
wave,  as  we  have  seen  above,  does  not  appear  at  reproducible  times  at  the 
reflecting  wall  in  a  shock  tube  experiment .  Te  began  our  study  of  the 
experiment  represented  by  Figure  3  by  calculating  induction  times  as  a 


TUBE  (°K) 
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Fig.  9  —  Resolution  test  of  the  calculated  temperature  as  a  function  of  position 
at  a  time  after  the  merging  of  reaction  wave  and  shock  front 
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Fig.  10  —  Calculated  number  density  of  OH  as  a  function  of  time  1  mm  from 
the  reflecting  wall  for  the  strong  ignition  case 
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function  of  temperature  for  the  conditions  behind  reflected  shocks  given 
in  Table  II.  By  analogy  with  the  strong  ignition  case,  in  the  ideal  ex¬ 
periment  and  calculation  the  reaction  wave  should  start  at  the  reflecting 
wall  about  3000  ysec  after  shock  reflection.  In  fact,  the  particular 
experiment  shown  in  Figure  3  indicates  that  the  reaction  wave  is  first 
observed  726  usee  after  reflection,  about  a  factor  of  four  too  soon.  We 
further  observe  that  the  change  in  induction  time  with  temperature  is  much 
greater  in  Figure  11  than  in  Figure  5. 

While  numerical  simulations  extending  past  3000  ysec  and  several 
meters  are  possible,  we  felt  that  we  could  learn  just  as  much  at  lower 
cost  by  increasing  the  reflected  shock  temperature  but  remaining  in  the 
same  general  vicinity  of  pressure  and  temperature  shown  in  Figure  11. 

Thus  we  simulated  a  case  for  which  the  temperature  and  pressure  behind 
the  reflected  shock  are  1000  K  and  3.72  atm  instead  of  972  K  and  3.4  atm 
for  the  same  H2:02:Ar/8:2:90  mixture.  The  calculated  induction  time  for 
this  mixture  is  M.550  ysec. 

Figure  12  shows  the  position  of  the  reflected  shock  front,  reaction 
wave,  and  the  transmitted  detonation  as  a  function  of  time  for  this  case. 
In  this  calculation,  the  reaction  wave  is  defined  as  that  place  where  the 
ratio  of  OH  to  Argon  becomes  significant  and  energy  has  just  begun  to  be 
released.  We  note  first  that  the  incident  shock  reflection  occurs  at 
VL00  ys  into  the  calculation.  The  first  indication  of  energy  release 
occurs  at  the  time  850  ys  which  is  750  ys  after  shock  reflection.  By 
1120  ys,  it  appears  that  a  detonation  has  formed  in  the  region  behind  the 
shock.  Figure  13  shows  profiles  of  temperature  as  a  function  of  position 
for  selected  times  during  the  period  from  820  to  1220  ys.  We  see  that  in 
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TEMPERATURE  (°K) 

Fig.  11  —  Calculations  of  the  chemical  induction  time  as  a  function  of 
temperature  for  the  weak  ignition  case 
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Fig.  12  —  Calculated  position  of  the  reflected  shock  front,  reaction  wave,  transmitted 
detonation,  and  contact  discontinuity  as  a  function  of  time  for  the  weak  ignition  case 
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fact  ignition  has  not  started  at  the  wall  in  the  calculation,  but  at  a 
location  about  3.5  cm  from  the  wall. 

A  close  look  at  the  calculation  indicates  that  there  are  effects  pre¬ 
sent  which  may  cause  the  results  to  deviate  from  ideality.  One  effect  is 
numerical  "clipping"  [20]  which  causes  the  temperatures  near  the  reflected 
wall  to  fluctuate  slightly  from  the  ideal  solution.  This  is  a  small  effect, 
causing  a  perturbation  of  a  few  degrees  (about  0.05%)  for  an  extent  of  about 
five  cells.  Such  a  temperature  perturbation  had  a  negligible  effect  on  the 
strong  ignition  problem  (Figure  5)  but  had  a  noticeable  effect  on  mild  ig¬ 
nition  (Figure  11).  Another  effect  observed  in  the  calculation  is  the 
presence  of  small  but  finite  velocities  behind  the  reflected  shock  during 
the  reflection  period.  As  shown  in  Figures  6  and  14,  these  velocities  start 

3 

out  at  fairly  large  amplitudes  (^5x10  cm/s)  and  then  die  out. 

Thus,  although  this  is  a  one-dimensional  calculation  and  we  therefore 
cannot  resolve  transverse  perturbations  and  boundary  layers,  we  find  that 
in  the  weak  ignition  case  the  system  is  extremely  sensitive  to  non-idealities 
and  small  perturbations.  We  have  thus  introduced  essentially  the  same  per¬ 
turbation  in  both  the  strong  and  the  weak  ignition  calculations.  In  the 
strong  ignition  case  we  have  obtained  quantitative  agreement  with  the  ex¬ 
periment  and  in  the  weak  ignition  case  we  obtain  agreement  qualitatively 
similar  to  that  observed  experimentally.  It  is  also  clear  that  the  ignition 
variability  in  the  weak  ignition  cases  is  adequately  reflected  in  the  much 

steeper  induction  time  versus  temperature  curve  of  Figure  11  relative  to 

3 

Figure  5.  In  the  weak  cases,  fluctuation  velocites  of  5x10  cm/sec  corre¬ 
spond  to  temperature  fluctuations  of  10°K  or  greater.  These  1%  fluctuations 
correspond  to  100%  variations  in  the  extremely  sensitive  chemical  induction 
time. 
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VI.  Discussion 

There  are  a  number  of  factors  that  might  help  explain  the  behavior 
observed  in  a  weak  ignition  experiments.  We  know  that  the  fluid  in 
the  shock  tube  behind  both  an  incident  and  reflected  shocks  is  not 
uniform.  Boundary  layers  and  transverse  waves  exist  which  perturb 
the  conditions  on  the  centerline.  Nonuniformities  in  the  incident 
shock  may  cause  it  to  hit  the  reflecting  wall  at  different  times  at 
each  position  in  the  wall.  Thus  shock  reflection  and  focusing 
effects  may  cause  uneven  ignition.  In  the  case  of  reflected  shocks, 
there  may  be  reflections  arising  from  the  interaction  of  the  reflect¬ 
ed  wave  and  the  contact  discontinuity  formed  at  the  onset  of  the  ex¬ 
periment.  Thus  when  the  reflected  wave  travels  long  times  and  dis¬ 
tances  before  the  reaction  wave  starts,  there  is  more  chance  for  the 
medium  to  be  perturbed.  This  is  especially  true  when,  as  we  see  in 
Figure  11,  the  medium  is  particularly  sensitive  to  perturbations. 

We  have  seen  that  there  are  deviations  from  ideality  in  the  one¬ 
dimensional  calculations  used  to  model  both  weak  and  strong  ignition. 
Large  oscillating  fluid  velocities  appear  near  the  wall  when  the  shock 
reflects.  Even  though  these  are  quickly  damped,  there  is  always  a  re¬ 
sidual,  nonzero  velocity  behind  the  reflected  shock.  There  is  also  a  de¬ 
viation  in  the  temperature  from  exactly  1000  K  at  the  wall,  the  tempera¬ 
ture  may  be  995  K  for  the  first  few  cells.  This  latter  effect  mimics 
the  effects  of  heat  loss  to  a  wall  due  to  thermal  conduction.  In  the 
case  of  mild  ignition,  where  the  induction  time  is  sensitive  to  pertur¬ 
bations,  we  expect  that  these  deviations  from  non-ideality  might  have  a 
noticeable  effect  on  our  answers. 


In  the  companion  paper  which  follows  this  we  have  analyzed  the  hydrogen- 
oxygen  system  with  respect  to  the  sensitivity  of  the  induction  time  to  sound 
wave  and  entropy  wave  perturbations.  This  had  led  to  the  development  of  a 


model  which  couples  the  chemical  kinetic  sensitivity  with  the  hydrodynamic 
fluctuations  to  quantify  the  system's  behavior.  This  model  clarifies  much 
of  the  weak  ignition  variability,  and  it  provides  a  scientific  basis  for  an 
improved  global  induction  time  model  [31,32]. 
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